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Abstract 

The search for a universal solution of the equations of motion for a 
satellite orbiting an oblate planet is a subject that has merited great 
interest because of its theoretical and practical implications. Here, a 
complete first-order perturbation solution, including the effects of the 
J 2 terms in the planet's potential, is given in terms of standard orbital 
parameters. The simple formulas provide a fast method for predict- 
ing satellite orbits that is more accurate than the two-bodv formulas. 
These predictions are shown to agree well with those of a completely 
numerical code and with actual satellite data. Also, in an appendix, it 
is rigorously proven that a satellite having negative mechanical energy 
remains for all time within a spherical annulus with radii approximately 
equal to the perigee and apogee of its initial osculating ellipse. 
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1 Introduction 



A characteristic feature of practical orbit prediction is that the engineer may deal with 
numerous satellites in a great variety of orbits. Under these circumstances analytical relations 
which can quickly approximate an orbit may be far superior to large numerical programs. 
While many analytical models have been developed for the artificial satellite age, most are 
not used in practical orbit prediction because they violate one or more of the following 
principles: 

• The method should provide a solution that is significantly more accurate than the 
two-body solution. 

• The real physical effects of the orbit should be easily distinguishable in the solution. 

• The solution should be universal; it should be valid for all orbital parameters. 

The problem of predicting the motion of a satellite perturbed only by the oblateness of the 
planet has received considerable attention following the first launchings of artificial satellites 
about the Earth. Some of the studies of this problem by means of general perturbation 
theories are listed at the end of this paper. Techniques have involved expansions in powers 
of \J~T2' averaging processes, the use of spheroidal coordinates, and the edifice of Hamiltonian 
mechanics. It is not the intention of this present paper to compare the various methodologies 
used. Suffice it to say that many researchers believe a solution which embodies all of the 
above principles was not achieved (e.g., see Taff). 

The basic procedure used in this paper to solve the differential equations of motion is 
the perturbation technique known as the Method of Strained Coordinates. This technique 
was first applied to the title problem by Brenner. Latta, and Weisfield. Using a mean orbital 
plane to specify an arbitrary orbit, they were only able to obtain a partial solution (e.g., the 
eccentricity was assumed small and initial conditions were not considered). 



9 






Here we use coordinates in the true orbital plane to cast the differential equations into a 
simplified form, as was originally done by Struble. 

2 Orbital Kinematics 

Figure 1 shows the usual reference system of spherical coordinates (r,a,0). The radial 
distance r is measured from the center of the planet 0 to the satellite S. The line O7 is in 
a direction fixed with respect to an inertial coordinate system. The right ascension a is the 
angle measured in the planet’s equatorial plane eastward from the line O7. The declination 
or latitude f3 is the angle measured northward from the equator. The position vector r of 
the satellite in the spherical coordinate system is 

r = r(cos q cos /3)b] + r(sin q cos /?)b 2 4- r(sin /3)b 3 (1 ) 

where (bi,b2,b 3 ) are orthonormal base vectors fixed in the directions shown. 

We can also locate the satellite by its polar coordinates (r. 0) within a (possibly rotating) 
orbital plane that instantaneously contains its position and velocity vectors. Here 6 is the 
argument of latitude, i.e., the angle measured in the orbital plane from the ascending node to 
the satellite. The orbital plane is inclined at an angle i to the equatorial plane and intersects 
the equatorial plane in the line of nodes, making an angle ft with the 0~j line. 

% 

We introduce another orthonormal set of base vectors (B].B2,B 3 ) which move with the 
satellite so that B, is in the direction of the position vector r. B 2 is also in the orbital plane, 
and B 3 = Bi x B 2 . The basis (bj , b 2 . b 3 ) may be transformed into the basis (B] . B 2 . B 3 ) by 
a succession of three rotations. First the basis (b],b 2 .b 3 ) is rotated about the b 3 direction 
by the angle ft, next the basis is rotated about the new 1-direction by the angle i, and 
finally the basis is again rotated about the new 3-direction by the angle 0. The two sets of 
base vectors are related by the product of the rotation matrices representing each successive 
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rotation (as explained in the book by Danielson): 



' B, ' 




cos 9 


sin 9 


0 ' 




' 1 


0 


0 * 




cosff 


sin 0 


0 ' 




b. 1 


b 2 


= 


— sin 9 


cos 9 


0 




0 


cos i 


sin i 




— sin Q 


cos f l 


0 




b 2 


b 3 




0 


0 


1 




0 


— sin i 


cos i 




0 


0 


1 




^3 



or 



B, ' 




cos 9 cos Q 


— sin 9 cos i sin D 


cos 9 sin Q. + sin 9 cos i cos fi 


sin 9 sin i 




' t>, ‘ 


b 2 


= 


— sin 9 cos 


— cos 9 cos i sin Q 


— sin 9 sin Q. + cos 9 cos i cos Q. 


cos 9 sin i 




b 2 


b 3 




sin i sin Q 


— sin i cos Q 


cos i 




^3 



The position vector r has only one component in the rotating basis: 



r = rB! 



(3) 



Using the first of equations (2). we obtain the components of r in the fixed basis: 

r = r(cos 9 cos fl — sin 9 cos i sin 0. )bj 

+ r(cos 9 sin Q 4- sin 9 cos i cos D)b 2 4- r(sin 6 sin t)b 3 (4) 

Equating the components of equations (1) and (4). we can obtain the following relations 
among the angles (a, 3) of the spherical coordinate system and the astronomical angles 

sin 3 = sin 6 sin i (5) 

cos 3 = cos 9 sec( o — V . ) 

The velocity dr/dt of the satellite is obtained by differentiating (3) with respect to the 
time t: 

dr dr r/R. 

( 6 ) 



dr dr c/B] 

Tt = li^ + r — 



Since the orbital plane must contain the velocity vector, we have to enforce 

d B, 



di 



B 3 = 0 



(") 



Substitution of equations (2) into equation (7) leads to a relationship which uncouples the 
equations for f 1(9) and i(9): 

dQ. tan 9 di 

(3) 



d9 sin i d9 
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The velocity (6) can then be written 



dr 

Tt 



dr d6 



' „ • di ' 

1 4- tan 9 cot i — 



B 



( 9 ) 



In the following part of this paper, we will obtain expressions for r(9), i{6), 0(0). and 
dt/d0(6). The position and velocity vectors of the satellite then may be calculated from 
the formulas in this section. The classical orbital elements p, e, and lj are the semilatus 
rectum, eccentricity, and argument of perigee of the instantaneous (osculating) conic section 
determined by the position and velocity vectors. If needed, p(0), e(0), and u >(8) can be 
obtained from our solution r(6) and dt/d6(6): 



P = 



(i) : 



e cos(# — u,’) = 1 






3 Equations of Motion 



The expressions for the kinetic and potential energies per unit mass of a satellite orbiting 
around an oblate planet are respectively. 

2 / , \ 2 / , \ 2 * 






d i) ^ 2 { d i) 



V = 



GM 






2r 2 



( 10 ) 



( 11 ) 



where G is the gravitational constant, M is the mass of the planet, R is the equatorial radius 
of the planet, and J 2 is the constant coefficient of the spherical harmonic of degree 2 and 
order 0 in the planet’s gravitational field. Substitution of these equations into Lagrange’s 
equations 

dd(T-V) „ 

T--T tt * U - M = o g = r,a. or p 

* a(S) 



0 



results in the following equations of motion: 



<Pr fd3V 2 0 ( da \ 2 dV 

777 ~ r 1 ^ 7 ) - rcos P\^7) =-fr 



( 12 ) 



d_ 

dt 




= 0 



d_ 

dt 




+ r 2 sin 3 cos 3 




2 



dV 

d3 



(13) 



Initial conditions are established by requiring that at the initial time to the orbital pa- 
rameters of the usual two-body orbit, the conic section determined by the initial position and 
velocity vectors, are known. The actual orbit is then tangent to this initial instantaneous 
conic section at t 0 (see Figure 1). Equating the initial position and velocity vectors given by 
equations (3) and (9) to the two-body expressions, we obtain 



r(t 0 ) — t 


Po 


(14) 


1 i- 6 0 COS((7 0 — uu q ) 




dr 


toho sin(#o — ^o) 




dJ Uo) ‘ 


Po 


(15) 


(M = — 


^0 

J "1 


(16) 


rl 


1 -f tan 0 O cot ?o^(0o)J 






O 

II 

O 


( 17 ) 


Wo) = 


(18) 



Here h 0 = \J 6 ' M po is the initial value of the satellite’s specific angular momentum about the 
center of the planet, and the subscript 0 on a symbol denotes that the parameter is evaluated 
at the initial time t 0 . 

We immediately have two integrals of the equations of motion: 



T + \ ’ = constant 



,da 



r 2 cos 2 3— = constant 

dt 



(19) 

( 20 ) 
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Equation (19) simply states that the mechanical energy of the satellite remains constant. 
Now, from equations (1) and (16) 



2 2 ^do dr , 

r cos p— = r x — • b 3 = h 0 cos x 0 
at at 



( 21 ) 



Equation (21) simply states that the component along the polar axis of the specific angular 
momentum of the satellite remains constant. Inserting equations (3) and (9) into equation 
(21), we obtain 



dt_ 

d6 



r 2 cos i 



' n . di ' 

1 + tan 0 cot x — 



( 22 ) 



/?0 COS Iq 

This allows the independent variable to be changed from t to 0. 

Letting u = po/r, and using equations (5). (21). and (22). we can rewrite the remaining 
equations of motion ( 1 2)— ( 13): 



di — 2Ju sin 0 cos 0 sin ? cos 2 f 
dO — — I- 2Ju sin 2 0 cos 3 i 

COS * 



d 2 u 

d(P 



+ u = 



cos 2 i J cos 2 i 



& 
d 2 u 



+ 



d 11 

u 2 ( 1 — 3 sin 2 0 sin 2 i) + 2u— sin 0 cos 0(1—3 cos 2 i) 

dO 



du' 



— 4 u sin 2 0 cos 2 i — 2 [ — ) sin 2 0 cos 2 i 
dO 2 l dO 



4 J 2 u sin 3 9 cos 6 i 



du t d f du\ 2 

“3r os<,(2 + sin ,) + Te{ u Te) sin6cos [ 

The terms in (24) with d 2 u/d0 2 can be combined, yielding the equivalent equation 

d 2 u ( cos 2 i J cos 2 i , ~ 

+ u = < — 1 — [u (1 + sin 0(7 cos x — 3)) 

a0 2 c 2 c l 



(24) 



du o ( du \ o o , 

-f 2u— sin 0 cos 0(1 — 3 cos 2 i) — 2 — sin 2 0 cos 2 1 ] 
du \ du J 



+ 



4J 2 u sin 3 0 cos 6 i 



2 • 2 . du 2 ( du\ 2 

u sin 0 cos x — u— cos 0(2 + sin x) — — sin 0 cos i 

dO \ dO I 



( 4 Ju sin 2 0 cos 4 x 4 J 2 u 2 sin 4 0 cos 8 ? 

v H ^ 1 : 



( 25 ) 



Here we have introduced the shorthand notation c = cosio, 5 = sin ?o. J = 12p^. 

4 Perturbation Procedure 



The differential equations (23 )— (24 ) are coupled by the nonlinear terms and apparently 
cannot be solved analytically. If we expand the right sides of (23) and (25) in a Taylor series 
expansion in powers of J and retain only terms up to order J 2 , the equations simplify to 

di — 2Jus\n0cos6s'micos 3 i 4 J 2 22 2 sin 2 cos ' 7 . 



d6 

d 2 u 

Jo 2 



+ 



cos 2 i J cos 2 i f —4 22 sin 2 6 cos 4 i 
+ u = — — + — 



sin 3 6 cos 6 + 0(J 3 ) 



+ u 2 [l + sin 2 0(1 cos 2 i — 3)] 



(26) 



+ s i n $ cos 0(1 — 3 cos 2 7 ) — 2 ( — ) sin 2 6 cos 2 ? 



dO 



-} 



(27) 



+ 



4 J 2 u sin 2 6 cos 6 2 



|u 2 [— 1 + 3 sin 2 0( 1—2 cos 2 i)] + 



3u sin 2 0 cos 4 i 



+ 22 — sin 6 cos 0[1 cos 2 i — 5] + ( - 37 ^ sin 2 0 cos 2 2 } + 0(J 3 ) 
do \ do I ) 



dO 



Here the term in the O symbols indicates that, for all sufficiently small J , the error is less 
than a constant times J 3 . The equations (26)— (27) are identical to those used as the starting 
point in the analysis of Eckstein, et al. 

It is reasonable to expect that the solution for u will be arbitrarily close to the two body 
solution. 1 + e o cos (0 — ^’ 0)1 when J is close to zero. This assumption is consistent with 
letting 

22 = 1 + to cos V 4" J + J 2 222 + . . . (28) 

y = 6 — u,'o + Jyi + J 2 t /2 + • • • (29) 

2 = 2 0 ~ b 1 + i/ 2 72 + • • • (30) 

An algorithm for the perturbation procedure is: 



Lei n = 1 

Substitute expressions (2b)-(30) into the equations of motion (26)-(27) 
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Equate the coefficients of J n 

Choose the arbitrary constants so secular terms will not arise. 

Solve for the n th order solution 
Satisfy the initial conditions (If)-(18) 

Iterate on n 

The calculations were carried out with the symbolic manipulation program MACSYMA. 
In this paper we only briefly outline these calculations; for more details see the theses of 
Sagovac and Snider. 

Beginning by substituting equations (28) and (30) into (26). and equating the terms 
multiplied by J, we obtain 

^ = —sc sin 26 - sin(y + 26) + ^psin (y - 26) (31) 

A solution to this equation is 

i\ = 77 cos 26 + cos (y + 20) + cos (y - 20) + 7\'i cos(2i / - 20) + A' 2 (32) 

2 6 2 

The last two terms may be added because they are to lowest order homogenous solutions 
to equation (30). The term multiplied by the constant I\\ was added to eliminate secular 
terms in z 2 ; note that differentiating this term with respect to 0 produces terms multiplied 
by J, from equation (29). The constant I\ 2 was added to satisfy the initial condition (IT), 
which implies that A(0 O ) = 0 so 

A' 2 = — 77 cos 20 o r— cos(30o — u/’o) — — — — cos(0 o + ^’ 0 ) — h\ cos 2u,' 0 

2 6 2 

Substituting equations (28)-(30) and (32) into (27). and equating terms multiplied by J 
yields 

(Pu, 3s 2 ■y ( bs 2 \ 1 r 9 9 

+ ui = 1 - — + t 2 0 ( — — + 1 J + -[(2 + 5e^)s 2 - 2 cq] cos 20 

e 2 £ lo£ 2 

+ (-9s 2 + 8) cos 2y + ^(1 Is 2 - 6) cos(y + 20) + —^(3 s 2 - 2) cos (2 y + 26) (33) 

4 o 24 
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+ 



- 2 ) - 

O c 



cos(2 y - 28) - + e 0 + 4 - 5s 2 'j cos y + sin y 



In the above equation, the cos y and siny terms would produce secular terms ^ sin j/ and 
Ocosy in U\. The choice dyi/d8 = 5s 2 /2 — 2 will eliminate these possibilities. Integrating 
yields 

yi = (t - 2 ) {e ~ 6o) + A ' 3[sin(2y - 2e) + sin 2u; ° ] (34) 

The term multiplied by A3 was added to eliminate secular terms in ti 2 . The constant terms 
in (34) were added to satisfy the initial condition y(9 0 ) = 8 0 — u 0 . 

A solution to lowest order of equation (33) is then 



3s 2 2 

u 1 = 1 — + e o 



+ ^(9s 2 -8)cos2i/ + ^(-ll5 2 + 6)cos(y + 20) + ^(-3s 2 +2)cos(2j/ + 2t?) (35) 

2s I\ j 



—os 



1 



+ 1 I + — [— s 2 (2 + 5eo) + 2 cq] cos 



+ 



f(3s 2 -2)- 
o c 



2s Ko 

cos(2 y — 26) 1 + J\ 4 cos (y — 26) 

c 



+ A' 5 cos (y - 6 Q + u.’ 0 ) + A' 6 sin(y - 8 0 + u; 0 ) 



The term multiplied by I\ 4 was added to eliminate secular terms in u 2 . The terms multiplied 
by A" 5 and Ag were added to satisfy the initial conditions (14)— (16). 

With all terms in place to deal with secular terms, the calculations are continued by 
substituting equations (28)-(30). (32). (34). and (35) into (26) and equating terms multiplied 
bv J 2 : 

sc£g( 1 5s 2 — 14) 



di2 

~d6 



A'i + 



sin(2y - 28) + . . . 



(36) 



24(5s 2 - 4) 

We have for brevity only indicated on the right side of equation (36) the term that would 
produce secular terms in ? 2 . Removal of this term by making its coefficient zero determines 
K\. Equation (36) is then integrated to determine ? 2 . 

Continuing the procedure by equating the terms multiplied by J 2 in the expansion of 
equation (27) determines y 2 . A3, and A' 4 . Final values of all the constants are listed in 
Appendix I. 
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Knowing the solution for z(0), we can determine Q(0) by integrating equation ($) and 
applying the initial condition (18). The angle 0, which increases continuously from an initial 
value 0 O , may be related to the time t by numerically integrating ( 22 ). 

5 Solution 



Here we assemble the complete solution: 

3s 2 



= po/|l + e 0 cos y + J 



1 7T + e o 1 - 



5s 2 ' 



+ — ( — (2 + 5cq)s 2 + 2 £q) cos 20 



2 2 

+ ~(9s 2 — 8 ) cos 2 y 4 * 7 -— ( — 1 Is 2 + 6 ) cos(y + 2 0 ) + ~("“3s 2 + 2 ) cos( 2 y + 26) 
12 24 24 



+ -^(3s 2 - 2 )cos( 2 y- 20 ) 
8 



+ 



e 0 [15(2 + ejj;)s 4 - 14(4 + e^s 2 + 24] sin [y- (5s 2 - 4)] sin[0 + u? 0 ] 

12(5s 2 — 4) 

e 2 s 2 (15s 2 - 14) sin [y (5s 2 -4)] sin 2 ^ 0 - y (5 s 2 - 4) 



6(5s 2 — 4) 



e 2 s 2 

tt - cos(y — 9o + 3^’o) 

lb 



2 2 2 

■f r~(3s 2 — 2) cos (y — 30 o + 3u; 0 ) yr cos(y — 50o + 3o.’ 0 ) 

24 lb 



+ — (3s 2 — 2) cos (y — 20q + 2a.’o) — 



3€qS 2 



cos (y — 40q + 2ll.’ 0 ) 



— ~r(s 2 + l)cos(y -I- 2 ^ 0 ) -I- -[(—2 -I- 5cq)s 2 — 2 cq] cos (y + 0 O + ^’ 0 ) 

4 O 



1 



1 



+ ”[(6 4- 5eo)s 2 — 4(1 + Cq )] cos (y — 6 0 + a-’o) 

+ — [—(14 + 5cq)s 2 + 2e 2 ] cos(y — 30 o + w 0 ) 

2 2 

+ ~(9s 2 — 4) cos(y + 30o — <^o) + ~ s? + 5) cos(y + 0o — <^o) 

4o o 



+ 7 t( — 5s 2 4 4) cos (y — 0o — ^’ 0 ) 
lb 



+ y (2s 2 - 1) cos(y + 20 o ) + j(-3s 2 + 1) cos(y - 20 o ) + + 2 ) cos V 



£0 , 

4 



4-eoS 2 cos(0q + u>o) 4 — ~ — cos(30q — <*-’o) + s 2 cos 20q] | + PqO(J^ , J 2 6) 



(37) 
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y — 6 — mo + J (~2 2^ (6 — 6 0 ) 



J4 j 


f (-75s 6 + 260s 4 - 296s 2 + 112) sin 


[f (5 s ! -4)] 


| COS | 


[2-.-o - f (5s 1 - 4)] 


' 24(5s 2 - 4)1 


1 

CN 

iO, 



+ J0s 2 ( — 15s 2 + 14)(15s 2 — 13)cos2u;o| 4- (15s 2 — 13) cos(0 o + w 0 ) 

2 2 
+ ^-(15s 2 - 13)cos(30 o - w 0 ) + ^-(15s 2 - 13)cos20 o 



+ ^:[5(9e 2 + 34 )s 4 + 4(9e 2 - 34)s 2 - 56e 2 ]} + 0(J\ J 3 6) (38) 



i = ? 0 4" scj i ~~ cos 20 ~~ cos(i/ -f* 20) 

[2 6 

e 0 eg(-15s 2 + 14) sin W (os 2 - 4)] sin [2a? 0 - y (5s 2 - 4)] 

+ ? co s („ - 20) + ' 12(5 f V- ~i) 



— — cos 26 o — — cos(30o — ~'o ) — ~ cos(0o + wo) 1 + 0( J 2 , J 3 6) 



(39) 



— fig + cj 



0o - 6 + ^ sin 26 - e 0 sin y + ^ sm(y + 26) - ^ s in(t/ - 26) - ^ sin 2 6 0 
2 o 2 1 



+c 0 sin(0 o - w 0 ) - y sin(30 o - w 0 ) - y sin(0 o + w 0 ) 



cJe 2 j 


|- 2(15s 4 — 45s 2 + 28) sin 


f (5s 2 -4)^ 


cos 


2w 0 - f(5^ 2 - 4) 


12(5s 2 - 4)1 


l (5s 2 -4) 



+J0s 2 (15s 2 — 14) cos 2^ 0 j -f cJ 2 6 — e 0 s 2 cos(0 o + u.’ 0 ) — cos(30 o — w 0 ) 



— s 2 cos 20 o + 24 C”' 5 ' 2 ~ ^ ^ ^ 



k+0(J 2 .J 3 6 ) 



(40) 



1 r8 f 

t — t Q + — / r 2 <l+J 
flQ J 6 q l 



(-3s 2 + 2) 



cos 2 6 + e 0 (s 2 - 1) 
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( 41 ) 



cos y + — ■ ^ cos (y + 20) + — — ~ - — — cos(y - 26) 

o 1 

Co5 2 (155 2 - 14) sin y (5s 2 — 4)] sin 2^ 0 - y (5s 2 — 4)] 

12(5s 2 — 4) 

2 2 2 

+5 2 — 1 + — cos 20q 4 cos (3^0 — <^'o) T — ~— cos(#o “I" <*?o) 

2 0 2 



Po. 



>d0+^O{J 2 ,J 2 0) 

ho 



In obtaining the equations (37)— (4 1 ), use has been made of trigonometric formulas 
to simplify terms containing the factor 5s 2 — 4 in the denominator. In the form given, 
these terms can clearly be seen to approach a finite limit at the “critical inclination'' 
to — sin' - 1 yfe/b = 63 0 26' or 116°34'. Hence the solution is actually valid for all values 
of io . If |?o — sin -1 ^/4/5| < J, the formulas ( 37 )— ( 4 1 ) can still be used by letting os 2 — 4 = J , 
or the limiting forms for i 0 — > sin -1 can be used. 

To check the solution, we can see if the specific mechanical energy (18) of the satellite 
remains constant. Substitution of the solution (36)— (37) into equation (10) plus (11) yields 



T + V = 



GM{ 1 — eg) GMJ 2 R 2 {\ -3sin 2 i 0 ) GM 



2p 0 



Or*3 



+ 



Po 



-0(J 2 ) 



The right side is easily recognized as the value of the specific mechanical energy at the initial 
time to. 

As a further check on the solution, we can see if it reduces to our previous results for 
equatorial and polar orbits, obtained by completely separate derivations (Danielson and 
Snider, 1989). Setting i 0 = 0 and using the independent variable a measured from the line 
O 7, we find that equations (37 )— (4 1 ) reduce to equations ( 1 8 )— ( 22 ) of our previous paper. 
Setting i 0 = 7r/2 and using the expansion cos (y + Jk) % cos y — Jksiny, we find that 
equations (37)— (41 ) reduce to equations (38)— (4 1 ) of our previous paper. 

Comparing the terms in the O-symbols, we see that the relative error in equation (41) 
may be greater than that of equations (37)-(40). Since the underlined terms in equations 
(37)-(40) are of this same order of magnitude, we can drop the underlined terms except 
when (37)— (38) are used to calculate r in equation (41). The relative error of our solution 
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will then still be of order (0 — &o)J 2 ■ 

If we retain only the two-body solution, the relative error terms will be of the order 
(0 — 6 0 )J. Here the error in our solution, as compared to the exact solution of the equations 
of motion, should be of the order J times the error in the two-body solution (for an Earth 
satellite J < .0015). 

6 Comparison of Perturbation, Two-Body, Numerical, and Mea- 
sured Solutions 



In this section we compare the preceding perturbation solution, the two-bodv solution, a 
completely numerical solution of the differentia] equations, and actual measured satellite 
data; for more comparisons see the thesis of Krambeck. The difference between the position 
vector r determined by the numerical integration code or measured data and the position 
vector r re f calculated from our perturbation solution or the two-body solution is the error 
Ar: 



Ar = r - r re f 



If the errors ( Ar. SO, Si, Aft) in the orbital parameters (r. 0, i , ft) are small, we can estimate 



Ar from equation (4) and the linear approximation 



dr dr dr dr 



(42) 



It is customary to decompose Ar into components (<5) , <5 2 . S3) along the moving triad (B] . B 2 , B 3 ): 



Ar — <5]B] •+■ 62 B 2 + 63 B 3 



The component 61 is called the radial error, 62 is the down track error, and S3 is the cross 
track error. Applying (42) to equation (4), and expressing the base vectors (bi,b 2 .b 3 ) in 
terms of (Bj,B 2 ,B 3 ), we obtain the following approximations: 

<5i%Ar. 62 % r(S 6 + cos i Aft) . S 3 % r(sin OSi — cos 8 sin i Aft) (43) 
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We obtained the numerical integration code UTOPIA from the Colorado Center for 
Astrodynamics Research located on the campus of the University of Colorado. The code 
was specialized to the differential equations used in this paper. We compared the solutions 
for an earth satellite with the following initial conditions: 

r 0 = 7,386.18 km 
Co = .003991 
0 O = 104.05° 
u.' 0 = 224.38° 
i 0 = 90.03° 
fio = 322.63° 

t o = 0 

These initial conditions represent an essentially polar orbit at an altitude of approximately 
1000 kilometers and period about l| hours. For this satellite the perturbation and numerical 
orbits match extremely well while the two-body orbit is grossly erroneous. The magnitude of 
the error in r is shown in Figure 2. Note that the relative error in our perturbation solution 
is 2.8 J 2 (0 — 0 o )i and that this error is 1.1 J times the error in the two-body solution. 

We obtained measured satellite data from the First Satellite Control Squadron located 
at Falcon Air Force Base. Colorado. A near earth satellite processed the following initial 
conditions: 



r 0 = 


7, 776.58 km 


Co — 


.0003071 


II 

O 


149.14° 


u?0 — 


9.57° 


*0 — 


98.81° 
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37.10° 



flo — 

t 0 = 0000Z 26 July 1990 

Again, the perturbation orbit is far superior to the two-body orbit. The radial, down track, 
and cross track errors (^1,62,63) are shown in Figure 3. Note that although the perturbation 
solution produces only a small improvement in the radial error, this error is negligible in 
comparison to the down track error. 

7 Conclusions 

Our solution embodies the principles outlined in the introduction. The relative error of our 
solution is of order (0 — 6 0 )J 2 , which is a factor of J times the relative error of the two-body 
solution; our solution loses its validity after an angular change (0 — 6 0 ) of order 1/J 2 , which 
is a factor of j longer than the interval of validity of the two-body solution. Secondly, our 
solution is in terms of classical orbital elements; no transformation to an alternative non- 
physical set of elements is required. Finally, our solution is free of singularities for all values 
of the initial orbital parameters, including elliptic, parabolic, and hyperbolic orbits. 

Our formulas should agree closely with satellite orbits whose dominant perturbation is 
the planet’s oblateness. Of course, the effects of higher-order terms in these expansions, 
higher-order terms in the planet's potential, and of other perturbation forces may also be 
important. The formulas will have to be amended to include these additional effects. 

APPENDIX I: Values of the Constants K\-Kq 



A’> = 



cse 2 0 { — 15s 2 + 14) 
24 (5s 2 -4) 



r- 5C on 5ce o , on x sct o ta , \ , cse£(15s 2 — 14) n 

A 2 = — — cos 20 o — cos(30 o — Wo) — cos (00 + w’o) H 94^2 _ 4) — cos 



A’, = 



£q(— 75s 6 + 260s 4 - 296s 2 + 112) 
48(5s 2 — 4) 2 
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r , _ [1 5 (Cq + 2)s 4 — 14(eo 4- 4)s 2 + 24] 

4 60 24 ( 5^ 2 - 4 ) 

2 2 

A 5 = ~( — 9«s 2 4 8) cos(2#o — 2u^o) 4* 777 ( 3s 2 — 2) cos(4$o — 2t^o) 

12 24 

2 

~(co^ 2 4~ A4) cos(6 0 + u?o) 4* ~(s 2 ~~ 2) cos(3$o — &o) 4* ~( — 3s 2 4 2) cos 2^o 

o o 

— 1^[5(2 — eo)-® 2 + 2 cq] cos 2^0 + p(15eo + 18 )s 2 — (cq + 1) 

2 2 

A 6 = +-~(65 2 — 5) sin (26$ — 2u?o) 4- “(—3s 2 + 1) sin(4#o ~ 2u^o) 
o 12 



“h — [^o( — 4* 1) + 2 A 4 ] sin(#o + ^’o) 4- — (3s 2 — 2) sin(# 0 — ^0) 



-f — ( — 7^ 2 4- 2) sin(30 o — ^0) 4 — ~(— s 2 4 1) sin 2wJq 4" 7 [~ ( 06q 4- 2)s 2 4- 2 Cq] sin 20q 
o 4 0 



APPENDIX II: Rigorous Bounds on the Orbit 



It follows from (10)— ( 12) that 



T + V-l($ V + :£-^ + £J^ (I -3» n >3) 



2 \dt J 2 dt 2 2 r 

This can be rewritten in the form 



4r 3 




dr 



from whence it follows that 

d_ 

dr 



= 4(T+V)r + 2GM + 



GMJoR 2 



(3 sin 2 /3 — 1 ) 




< 4(7 + V)r + 2 GM + 2GM ^ R2 

r* 2 



Integrating from r(f 0 ) to r(t) yields 




r* | — I < 2(7 + V)r 2 + 2GM r - 



2 GMJ 2 R 2 



-h 2 0 + 



ZhjJ 2 R 2 

Po r o 



It follows that 



0 < 2{T + V)r 2 + 2GMr - h 2 {l - 



3 J 2 /? 2 
Poc 0 



(44) 



17 



When T + V < 0, the quadratic polynomial on the right side of (44) has the roots (exact 
values can be found from the quadratic formula) 



Tmin = T^-[l + 0(J 2 )} , r max = -£2_[1 + 0(J 2 ) } 

Hence a satellite having negative mechanical energy remains for all time within the spherical 
annulus r^n < r < r^x. Since the position vector is bounded, we can invoke the recurrence 
theorem; i.e., the satellite will come as close as desired to its initial position in a sufficiently 
long period of time (as shown by Poincare). Furthermore, we are guaranteed of the validity 
of supressing secular terms to describe the orbit via perturbation analysis. 
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Figure 1: Orbital geometry. 



Figure 2: Comparison of perturbation, two-body, and numerical orbits. 



Figure 3: Comparison of perturbation, two-body, and measured orbits. 
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